/*
 * CarbonBasedAblationMaterial.C
 *
 *  Created on: Dec 24, 2019
 *      Author: liuxiao
 *
 */

#include "CarbonBasedAblationMaterial.h"
#include<iostream>
#include "libmesh/quadrature.h"
#include "libmesh/system.h"
#include "libmesh/radial_basis_interpolation.h"
using std::cout;
using std::endl;
registerMooseObject("TrilobitaApp", CarbonBasedAblationMaterial);
template<>
InputParameters validParams<CarbonBasedAblationMaterial>()
{
	InputParameters params = validParams<Material>();
	params.addParam<Real>("qc", 0, "cold wall heat flux at Tw=0K");
	params.addParam<Real>("hr", 0, "enthalpy of recovery");
	params.addParam<Real>("Pe", 1, "pressure out of boundary layer,atm");
	params.addParam<Real>("Pe_O2", 0.21, "pressure of O2 out of boundary layer,1");
	params.addParam<Real>("Pe_N2", 0.79, "pressure of N2 out of boundary layer,1");
	params.addParam<Real>("M_e",0.29, " average molar mass of gas out of boundary layer,kg/mol");
	params.addParam<Real>("M_p",0.44, " average molar mass of pyrolysis gas,kg/mol");
	params.addParam<Real>("mp",0, " mass flux of pyrolysis gas,kg/m**2/s");
	params.addRequiredCoupledVar("Tw", "Coupled wall temperature");
	params.addParam<int>("gamma",0, "to choose laminar 0 or turbulence 1");
	params.addParam<int>("scala_model",1, "to choose scala model,0=slow reaction,1=fast reaction");

	return params;
}


CarbonBasedAblationMaterial::CarbonBasedAblationMaterial(const InputParameters & parameters) :
    	    						Material(parameters),
									_Pw_C1(declareProperty<Real>("Pw_C1")),
									_Pw_C2(declareProperty<Real>("Pw_C2")),
									_Pw_C3(declareProperty<Real>("Pw_C3")),
									_Pw_C4(declareProperty<Real>("Pw_C4")),
									_Pw_C5(declareProperty<Real>("Pw_C5")),
									_Pw_O2(declareProperty<Real>("Pw_O2")),
									_Pw_O(declareProperty<Real>("Pw_O")),
									_Pw_CO(declareProperty<Real>("Pw_CO")),
									_Pw_CO2(declareProperty<Real>("Pw_CO2")),
									_Pw_N2(declareProperty<Real>("Pw_N2")),
									_Pw_N(declareProperty<Real>("Pw_N")),
									_Pw_CN(declareProperty<Real>("Pw_CN")),
									_Pw_C2N(declareProperty<Real>("Pw_C2N")),
									_Cw_C1(declareProperty<Real>("Cw_C1")),
									_Cw_C2(declareProperty<Real>("Cw_C2")),
									_Cw_C3(declareProperty<Real>("Cw_C3")),
									_Cw_C4(declareProperty<Real>("Cw_C4")),
									_Cw_C5(declareProperty<Real>("Cw_C5")),
									_Cw_O2(declareProperty<Real>("Cw_O2")),
									_Cw_O(declareProperty<Real>("Cw_O")),
									_Cw_CO(declareProperty<Real>("Cw_CO")),
									_Cw_CO2(declareProperty<Real>("Cw_CO2")),
									_Cw_N2(declareProperty<Real>("Cw_N2")),
									_Cw_N(declareProperty<Real>("Cw_N")),
									_Cw_CN(declareProperty<Real>("Cw_CN")),
									_Cw_C2N(declareProperty<Real>("Cw_C2N")),
									_M_bar(declareProperty<Real>("M_bar")),
									_Bc(declareProperty<Real>("Bc")),
									_Bp(declareProperty<Real>("Bp")),
									_Bw(declareProperty<Real>("Bw")),
									_mc(declareProperty<Real>("mc")),
									_injection_coff(declareProperty<Real>("injection_coff")),
									_qc(getParam<Real>("qc")),
									_hr(getParam<Real>("hr")),
									_Pe(getParam<Real>("Pe")),
									_Pe_O2(getParam<Real>("Pe_O2")),
									_Pe_N2(getParam<Real>("Pe_N2")),
									 _M_e(getParam<Real>("M_e")),
									 _M_p(getParam<Real>("M_p")),
									 _mp(getParam<Real>("mp")),
									_Tw(coupledValue("Tw")),
									_gamma(getParam<int>("gamma")),
									_scala_model(getParam<int>("scala_model")),
									_Kp_C1(0),
									_Kp_C2(0),
									_Kp_C3(0),
									_Kp_C4(0),
									_Kp_C5(0),
									_Kp_O(0),
									_Kp_CO2(0),
									_Kp_CN(0),
									_Kp_C2N(0),
									_Kp_N(0),
									_M_C(0.012),
									_M_C1(0.012),
									_M_C2(0.024),
									_M_C3(0.036),
									_M_C4(0.048),
									_M_C5(0.060),
									_M_O2(0.032),
									_M_O(0.016),
									_M_CO(0.028),
									_M_CO2(0.044),
									_M_N2(0.028),
									_M_N(0.014),
									_M_CN(0.026),
									_M_C2N(0.046),
									_R(8.314)
{

}

void CarbonBasedAblationMaterial::compute_Kp(Real T)
{
	_Kp_C1=pow(10,(8.12-3.7219e4/T));
	_Kp_C2=pow(10,(9.64-4.2713e4/T));
	_Kp_C3=pow(10,(9.78-4.0481e4/T));
	_Kp_C4=pow(10,(10.013-4.889e4/T));
	_Kp_C5=pow(10,(10.464-4.927e4/T));
	_Kp_O=pow(10,(3.518-13392/T));
	_Kp_CO2=10*pow(10,(-4.482+1.4608e4/T));
	_Kp_CN=pow(10,(5.029-22288/T));
	_Kp_C2N=pow(10,(6.572-28450/T));
	_Kp_N=pow(10,(3.542-25176/T));
	}

Real CarbonBasedAblationMaterial::compute_Ci_from_P(Real P,Real Mi,Real M_bar)
{
	Real C_i;
	C_i=(P/_Pe)*(Mi/M_bar);
	return C_i;
	}

Real CarbonBasedAblationMaterial::compute_P_Ci(int i,Real Bw,Real phi,Real M_bar,Real T)
{
	std::vector<Real> P_Ci;
    std::vector<Real> beta;
    std::vector<Real> alpha;
    std::vector<Real> P_Ceq;
    std::vector<Real> M;
    Real Phol;
    Real Ahol;
    beta.push_back(0);
    beta.push_back(0.24);
    beta.push_back(0.5);
    beta.push_back(0.023);
    beta.push_back(0.25);
    beta.push_back(0.0019);
    alpha.push_back(0);
    alpha.push_back(beta[1]*M_bar*_Pe/sqrt(2*pi*_R*_M_C1*T));
    alpha.push_back(beta[2]*M_bar*_Pe/sqrt(2*pi*_R*_M_C2*T));
    alpha.push_back(beta[3]*M_bar*_Pe/sqrt(2*pi*_R*_M_C3*T));
    alpha.push_back(beta[4]*M_bar*_Pe/sqrt(2*pi*_R*_M_C4*T));
    alpha.push_back(beta[5]*M_bar*_Pe/sqrt(2*pi*_R*_M_C5*T));
    P_Ceq.push_back(0);
    P_Ceq.push_back(_Kp_C1);
    P_Ceq.push_back(_Kp_C2);
    P_Ceq.push_back(_Kp_C3);
    P_Ceq.push_back(_Kp_C4);
    P_Ceq.push_back(_Kp_C5);
    M.push_back(0);
    M.push_back(_M_C1);
    M.push_back(_M_C2);
    M.push_back(_M_C3);
    M.push_back(_M_C4);
    M.push_back(_M_C5);
    Phol=phi*_qc/_hr;
    Ahol=M_bar*_Pe;
    for (unsigned int i(0);i < 6; ++i)
    {
    	if(i==0)
    	{
    		P_Ci.push_back(0);
    	}
    	else
    	{
    		P_Ci.push_back(P_Ceq[i]/(1+(1+Bw)*Phol/(alpha[i]*Ahol/M[i])));
    	}

    }




    return P_Ci[i];

}

Real CarbonBasedAblationMaterial::compute_Pw_O2(Real Pw_CO,Real Pw_CO2,Real Pw_O,Real Bw,Real phi,Real M_bar,Real T)
{
	Real A1;
	Real B1;
	Real alpha1;
	Real A1_p;
	Real B1_p;
	Real alpha1_p;
	Real Phol;
	Real a;
	Real b;
	Real c;
	Real X;
	Real P_CO_CO2;
	Real P_O_CO2;
	Real alch;
	int f;
	if(_scala_model==0)
	{
		A1=2.18e5;
		B1=21558;
		alpha1=A1*exp(-B1/T);
	}
	else if(_scala_model==1)
	{
		A1=3.29e9;
		B1=22144;
		alpha1=A1*exp(-B1/T);
	}
	else
	{
		 mooseError("wrong scala_model value,must be 0 or 1");
	}
   A1_p=4.88e10;
   B1_p=42500;
   alpha1_p=A1_p*exp(-B1_p/T);
//   if(T<=3300)
//   {
	   f=0;
//   }
//   else
//   {
//	   f=1;
//   }

   if(abs(Pw_CO2<1e-10))
   {
	   Pw_CO2=1e-10;
   }
   P_CO_CO2=Pw_CO/Pw_CO2;
   Pw_O = Pw_O/Pw_CO2;
   alch=(2+P_CO_CO2+Pw_O)*_M_O2/(2+2*P_CO_CO2)/_M_C;
   Phol=phi*_qc/_hr;
   a=Phol*_M_O2*(1+Bw)/M_bar+f*alpha1_p*alch;
   b=(1-f)*alpha1*alch;
   c=-Phol*_M_O2*_Pe_O2/M_bar;
   if((pow(b,2)-4*a*c)>0&&abs(a)>1e-32)
   {
   X=(-b+sqrt(pow(b,2)-4*a*c))/(2*a);
   }
   else
   {
	   X=0;
	   mooseError("Pw_O2 has no solution!");
   }
   Real mc1=(1-f)*alpha1*X+f*alpha1_p*pow(X,2);

   return pow(X,2);

}
Real CarbonBasedAblationMaterial::compute_Pw_O(Real Pw_O2)
{
	Real Pw_O;
	Pw_O=Pw_O2*_Kp_O;
	return Pw_O;
}

Real CarbonBasedAblationMaterial::compute_Pw_CO(Real Pw_O2,Real Pw_O,Real M_bar,Real Bw)
{
	Real Ce_O;
	Real Pw_CO;
	Ce_O=compute_Ci_from_P(_Pe_O2,_M_O2,_M_e);
	Pw_CO=(Ce_O*M_bar*_Pe/(1+Bw)/_M_O-Pw_O-2*Pw_O2)/(1+2*sqrt(Pw_O2)*_Kp_CO2);
	return Pw_CO;
}

Real CarbonBasedAblationMaterial::compute_Pw_CO2(Real Pw_CO,Real Pw_O2)
{
	Real Pw_CO2;
	Pw_CO2=_Kp_CO2*Pw_CO*sqrt(Pw_O2);
	return Pw_CO2;
}

Real CarbonBasedAblationMaterial::compute_Pw_N2(Real Bw,Real phi,Real M_bar)
{
	Real Ce_N;
    Real a;
    Real b;
    Real c;
    Real Y;
    Ce_N=compute_Ci_from_P(_Pe_N2,_M_N2,_M_e);
    a=2;
    b=_Kp_N+_Kp_CN+_Kp_C2N;
    c=-_Pe*M_bar*Ce_N/_M_N/(1+Bw);
    if((pow(b,2)-4*a*c)>0)
      {
      Y=(-b+sqrt(pow(b,2)-4*a*c))/(2*a);

      }
      else
      {
   	   Y=0;
   	   mooseError("Pw_O2 has no solution!");
      }
      return pow(Y,2);
}
Real CarbonBasedAblationMaterial::compute_Pw_N(Real Pw_N2)
{
	Real Pw_N;
	Pw_N=_Kp_N*sqrt(Pw_N2);
	return Pw_N;
}

Real CarbonBasedAblationMaterial::compute_Pw_CN(Real Pw_N2)
{
	Real Pw_CN;
	Pw_CN=_Kp_CN*sqrt(Pw_N2);
	return Pw_CN;
}
Real CarbonBasedAblationMaterial::compute_Pw_C2N(Real Pw_N2)
{
	Real Pw_C2N;
	Pw_C2N=_Kp_C2N*sqrt(Pw_N2);
	return Pw_C2N;
}

Real CarbonBasedAblationMaterial::compute_M_bar(Real Cw_C1,Real Cw_C2,Real Cw_C3,Real Cw_C4,Real Cw_C5,Real Cw_O2,Real Cw_O,Real Cw_CO,Real Cw_CO2,Real Cw_N2,Real Cw_N,Real Cw_CN,Real Cw_C2N)
{
	Real M_bar;
	M_bar=1/(Cw_C1/_M_C1+Cw_C2/_M_C2+Cw_C3/_M_C3+Cw_C4/_M_C4+Cw_C5/_M_C5\
			+Cw_O/_M_O+Cw_O2/_M_O2+Cw_CO/_M_CO+Cw_CO2/_M_CO2+Cw_N/_M_N+Cw_N2/_M_N2\
			+Cw_CN/_M_CN+Cw_C2N/_M_C2N);
	return M_bar;

}

Real CarbonBasedAblationMaterial:: compute_carbon_w(Real Cw_C1,Real Cw_C2,Real Cw_C3,Real Cw_C4,Real Cw_C5,\
		                                            Real Cw_CO,Real Cw_CO2,Real Cw_CN,Real Cw_C2N)
{
	Real carbon_w;
	carbon_w=Cw_C1+Cw_C2+Cw_C3+Cw_C4+Cw_C5+Cw_CO*_M_C/_M_CO+Cw_CO2*_M_C/_M_CO2\
			 +Cw_CN*_M_C/_M_CN+Cw_C2N*_M_C/_M_C2N;
	return carbon_w;
}

Real CarbonBasedAblationMaterial::compute_Bc(Real carbon_w,Real Bp)
{

	Real Bc;
	Bc = (1+Bp)*carbon_w/(1-carbon_w);
	return Bc;
}

Real CarbonBasedAblationMaterial::compute_Bp(Real mp,Real phi)
{

	Real Bp;
	Bp = mp/(phi*_qc/_hr);
	return Bp;
}
Real CarbonBasedAblationMaterial::compute_mc(Real Bc,Real phi)
{
	 Real mc;
	 mc = Bc*phi*_qc/_hr;
	 return mc;
}


Real CarbonBasedAblationMaterial::compute_phi(Real Bc,Real M_bar,Real mc,Real mp,int gamma,Real T)
{
Real kesi= (pow((M_bar / _M_C), 0.26) * mc + pow((M_bar / _M_p), 0.26) * mp) *_hr/_qc;
Real phi;
if(T<3300)
{
	if(kesi<2.25)
	{
		phi=1-0.724*kesi+0.13*pow(kesi,2);
	}
	else
	{
		phi=0.04;
	}
}
else
{
     phi=1-(1-gamma)*0.58*(mc+mp)*_hr/_qc-gamma*2*(mc+mp)*_hr/_qc;
}

return phi;

}

void CarbonBasedAblationMaterial::iteration(Real T,int qp)
{
  Real Cw_C1_it=0;
  Real Cw_C2_it=0;
  Real Cw_C3_it=0;
  Real Cw_C4_it=0;
  Real Cw_C5_it=0;
  Real Cw_O2_it=0;
  Real Cw_O_it=compute_Ci_from_P(_Pe_O2,_M_O2,_M_e);
  Real Cw_CO_it=0;
  Real Cw_CO2_it=0;
  Real Cw_N2_it=compute_Ci_from_P(_Pe_N2,_M_N2,_M_e);
  Real Cw_N_it=0;
  Real Cw_CN_it=0;
  Real Cw_C2N_it=0;
  Real Pw_C1_it=0;
  Real Pw_C2_it=0;
  Real Pw_C3_it=0;
  Real Pw_C4_it=0;
  Real Pw_C5_it=0;
  Real Pw_O2_it=_Pe_O2;
  Real Pw_O_it=0;
  Real Pw_CO_it=0;
  Real Pw_CO2_it=0;
  Real Pw_N2_it=_Pe_N2;
  Real Pw_N_it=0;
  Real Pw_CN_it=0;
  Real Pw_C2N_it=0;
  Real carbon_w_it=0;
  Real mc_it=0;
  Real Bc_it0=0;
  Real Bc_it1=0;
  Real Bp_it=0;
  Real delta_Bc_it=0.0001;
  Real M_bar_it=_M_e;
  Real phi_it=1;
  Real Res=1e-6;
  Real Bw_it=0;
  int init=1;
  compute_Kp(T);
  while(init==1||abs(Bc_it1-Bc_it0)>Res)
  {
	  init=init+1;
	  delta_Bc_it=(Bc_it1-Bc_it0)/2;
	  Bc_it0=Bc_it0+delta_Bc_it;
	  Bp_it=compute_Bp(_mp,phi_it);
	  Bw_it=Bc_it0+Bp_it;
	  Pw_O2_it=compute_Pw_O2(Pw_CO_it,Pw_CO2_it,Pw_O_it,Bw_it,phi_it,M_bar_it,T);
	  Pw_O_it=compute_Pw_O(Pw_O2_it);
	  Pw_CO_it=compute_Pw_CO(Pw_O2_it,Pw_O_it,M_bar_it, Bw_it);
	  Pw_CO2_it=compute_Pw_CO2(Pw_CO_it,Pw_O2_it);
	  Cw_O2_it=compute_Ci_from_P(Pw_O2_it,_M_O2,M_bar_it);
	  Cw_O_it=compute_Ci_from_P(Pw_O_it,_M_O,M_bar_it);
	  Cw_CO_it=compute_Ci_from_P(Pw_CO_it,_M_CO,M_bar_it);
	  Cw_CO2_it=compute_Ci_from_P(Pw_CO2_it,_M_CO2,M_bar_it);
	  if(T>=3000)
	  {
		  Pw_C1_it=compute_P_Ci(1, Bw_it,phi_it,M_bar_it,T);
		  Pw_C2_it=compute_P_Ci(2, Bw_it,phi_it,M_bar_it,T);
		  Pw_C3_it=compute_P_Ci(3, Bw_it,phi_it,M_bar_it,T);
		  Pw_C4_it=compute_P_Ci(4, Bw_it,phi_it,M_bar_it,T);
		  Pw_C5_it=compute_P_Ci(5, Bw_it,phi_it,M_bar_it,T);
		  Pw_N2_it=compute_Pw_N2(Bw_it,phi_it,M_bar_it);
		  Pw_N_it=compute_Pw_N(Pw_N2_it);
		  Pw_CN_it=compute_Pw_CN(Pw_N2_it);
		  Pw_C2N_it=compute_Pw_C2N(Pw_N2_it);
		  Cw_C1_it=compute_Ci_from_P(Pw_C1_it,_M_C1,M_bar_it);
		  Cw_C2_it=compute_Ci_from_P(Pw_C2_it,_M_C2,M_bar_it);
		  Cw_C3_it=compute_Ci_from_P(Pw_C3_it,_M_C3,M_bar_it);
		  Cw_C4_it=compute_Ci_from_P(Pw_C4_it,_M_C4,M_bar_it);
		  Cw_C5_it=compute_Ci_from_P(Pw_C5_it,_M_C5,M_bar_it);
		  Cw_N2_it=compute_Ci_from_P(Pw_N2_it,_M_N2,M_bar_it);
		  Cw_N_it=compute_Ci_from_P(Pw_N_it,_M_N,M_bar_it);
		  Cw_CN_it=compute_Ci_from_P(Pw_CN_it,_M_CN,M_bar_it);
		  Cw_C2N_it=compute_Ci_from_P(Pw_C2N_it,_M_C2N,M_bar_it);
	  }
	  carbon_w_it=compute_carbon_w(Cw_C1_it,Cw_C2_it, Cw_C3_it, Cw_C4_it,Cw_C5_it,Cw_CO_it,Cw_CO2_it,Cw_CN_it,Cw_C2N_it);
	  Bc_it1=compute_Bc(carbon_w_it,Bp_it);
	  mc_it=compute_mc( Bc_it1, phi_it);
	  M_bar_it=compute_M_bar(Cw_C1_it,Cw_C2_it,Cw_C3_it,Cw_C4_it,Cw_C5_it,Cw_O2_it,Cw_O_it,Cw_CO_it,Cw_CO2_it,Cw_N2_it,Cw_N_it,Cw_CN_it,Cw_C2N_it);
	  phi_it=compute_phi(Bc_it1,M_bar_it,mc_it,_mp,_gamma,T);
	  Real P_CO_CO2=Pw_CO_it/(Pw_CO2_it+1e-10);
  }
//Cw_C1=Cw_C1_it;
//Cw_C2=Cw_C2_it;
//Cw_C3=Cw_C3_it;
//Cw_C4=Cw_C4_it;
//Cw_C5=Cw_C5_it;
//Cw_O2=Cw_O2_it;
//Cw_O=Cw_O_it;
//Cw_CO=Cw_CO_it;
//Cw_CO2=Cw_CO2_it;
//Cw_N2=Cw_N2_it;
//Cw_N=Cw_N_it;
//Cw_CN=Cw_CN_it;
//Cw_C2N=Cw_C2N_it;
//Bc=Bc_it0;
//std::cout<<"Bc="<<Bc<<std::endl;
//M_bar=M_bar_it;
//phi=phi_it;
_Cw_C1[qp]=Cw_C1_it;
_Cw_C1[qp]=Cw_C1_it;
_Cw_C2[qp]=Cw_C2_it;
_Cw_C3[qp]=Cw_C3_it;
_Cw_C4[qp]=Cw_C4_it;
_Cw_C5[qp]=Cw_C5_it;
_Cw_O2[qp]=Cw_O2_it;
_Cw_O[qp]=Cw_O_it;
_Cw_CO[qp]=Cw_CO_it;
_Cw_CO2[qp]=Cw_CO2_it;
_Cw_N2[qp]=Cw_N2_it;
_Cw_N[qp]=Cw_N_it;
_Cw_CN[qp]=Cw_CN_it;
_Cw_C2N[qp]=Cw_C2N_it;
_Bc[qp]=Bc_it0;
//std::cout<<"Bc="<<_Bc[qp]<<std::endl;
_M_bar[qp]=M_bar_it;
_injection_coff[qp]=phi_it;
}


void CarbonBasedAblationMaterial::computeProperties()
{


	for (unsigned int qp(0);qp < _qrule->n_points(); ++qp)
	{
		iteration(_Tw[qp],qp);
	}

}





